All Morphology: Wing Time Series

Introduction

When to use a time series vs. a linear regression?

If you have continuous target variable, then it is a regression problem. For instance, in flight trials we have the flight distance to predict, which is continuous. Hence this becomes a regression problem. About time series, when the datapoints are time dependent, then it becomes a time series problem. Each data point has an order and is, typically, related to the data points before and after by some underlying process.

In turn, a times series regression is a statistical method for predicting a future response based on the response history (known as autoregressive dynamics) and the transfer of dynamics from relevant predictors. Time series regression is commonly used for modeling and forecasting of economic, financial, and biological systems.

There are three concepts then to keep in mind:

• Stationary: if any trends or patterns in the time series are independent with time. (A linear regression equivalent to a time series at stationarity). There is strong and weak stationary.

• Non-Stationary: trends dependent in time in the mean or statistical deviation of the data points.

• Autocorrelation: “memory”; the degree to which time series values in period(t) are related to time series values in periods (t+1, t+2, t+3…). For this, we can test to see how long certain values last in a system.

There are also three words and their definitions to remember:

• variance = amplitude of the data points
• mean = the central value of the data points
• covariance = frequency of the data points

Questions about time series often asked include,

• Is there any trend, seasonality, or any outliers?
• Is there a long-run cycle or period unrelated to seasonality factors?
• Is there a constant variance over time?
• Are there any abrupt changes to either the level of the series or the variance?

So how can you answer these questions through a time series regression and code? And where do you begin to test stationarity or visualize autocorrelation? And finally, how do you handle non-stationarity?

1. First, begin by cleaning the data.

Cleaning Data and Sourcing Scripts

source_path = "~/Desktop/git_repositories/SBB-dispersal/avbernat_working_on/Rsrc/"

script_names = c("compare_models.R","regression_output.R", "clean_morph_data.R",
                 "AICprobabilities.R")

for (script in script_names) { 
  path = paste0(source_path, script)
  source(path) 
}

Read the data

2. Extract date information in order to prep it for converting it into a datetime object read by the xts() function.

3. Merge close dates together in order to have as equally spaced datapoints as possible. The acf() assumes that the data are regular spaced. The acf computes estimates of the autocovariance or autocorrelation function. We need this because it will tell us, if there are patterns, what the mathematical function of those patterns are.

data_list <- read_morph_data("data/allmorphology9.21.20.csv")
raw_data = data_list[[1]]
all_bugs = nrow(raw_data)
# data_long = data_list[[2]] # need to fix this 

# Remove individuals with torn wings first
raw_data$drop <- FALSE
for(row in 1:nrow(raw_data)){
    if(length(unlist(strsplit(strsplit(paste("test ", raw_data$notes[row], " test", sep=""), "torn")[[1]], "wing")))>2){
         #browser() 
         raw_data$drop[row] <- TRUE
         }
}
raw_data <- raw_data[raw_data$drop==FALSE,]
clean_bugs = nrow(raw_data)

cat("number of bugs with torn wings:", all_bugs - clean_bugs, "\n\n")

## number of bugs with torn wings: 141

# Datetime
raw_data$date <- paste(raw_data$month, raw_data$year, sep="/")
raw_data$datetime <- as.yearmon(raw_data$date, "%B/%Y")
raw_data$datetime <- as.factor(raw_data$datetime)
n_missing_dates = nrow(raw_data[is.na(raw_data$datetime),])

# merge May 2015 with April 2015 because very few bugs were collected in May 2015.
# then merge April 2013 with May 2013 to make the time datapoints more evenly distanced
raw_data$date[raw_data$date == "May/2015"] = "April/2015"
raw_data$date[raw_data$date == "May/2013"] = "April/2013"

# convert to yearmon object and then factor
raw_data$datetime <- as.yearmon(raw_data$date, "%B/%Y")
raw_data$datetime <- as.factor(raw_data$datetime)

cat("number of missing dates:", n_missing_dates, "\n\n")

## number of missing dates: 296

unique(raw_data$datetime)

##  [1] Apr 2013 Apr 2014 Apr 2015 Dec 2013 Aug 2017 Dec 2016 Sep 2018 May 2019
##  [9] Oct 2019 Feb 2020 <NA>    
## 10 Levels: Apr 2013 Dec 2013 Apr 2014 Apr 2015 Dec 2016 Aug 2017 ... Feb 2020

Looks like we have a measurement at least once per year, which is good.

Time Series: All Data

Wing Length

4. Remove missing wing and datetime values.

5. Compute the wing length averages and generate datetime objects using as.Date(). In order to be processed by the as.Date() function, all time objects need a date, month, and year.

6. Use the xts & zoo R Libraries to read cleanly index the data by a formal time object (collection_time).

7. Optional: include major events that occurred in the 2010’s. This could be events that you find biologically significant to your questions. E.g. how did major hurricanes in Florida impact average soapberry bug wing length across Florida?

Hurricane Source: https://en.wikipedia.org/wiki/List_of_Florida_hurricanes#2000–present

# remove NA dates
d = raw_data %>%
  filter(!is.na(wing), !is.na(datetime))

# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

# events
FL_major_hurr = c("Sep 2017 10", "Oct 2018 10") # Irma (s), Michael (n)
hurr_dates = as.Date(FL_major_hurr, "%b %Y %d")

events <- xts(c("Irma (S. FL)", "Michael (N. Fl)"), 
              hurr_dates)

8. Create xts-zoo object and plot the time series data. Add any events to the plot with addEventLines().

9. Calculate the Augmented Dickey-Fuller Test to test for stationarity where the null hypothesis is non-stationarity and the alternative hypothesis is stationarity.

10. Finally, use an ACF (autocorrelation function) and PCF (partial-autocorrelation function) plot to identify temporal dependence in the data. Autocorrelation measures the linear relaionship between lagged values of a time series.

\(R_s = Corr(x_t, x_{t+s})\) for lag \(s\).

There are two horizonal, blue, dashed lines as well. Those represent the significance threshold, where only the spikes that exceed this dashed line are considered significant. In other words, these plots describe how well the present value of the series is related with its past values. Since, none do significantly, then there is no AR (a present value of the time series cannot be obtained using previous values of the same time series).

check_stationarity = function(dx) {
  #  dx is your xts zoo time series object
  plot(dx)
  #addEventLines(events, pos=2, srt=90, col="red")
  
  par(mfrow=c(1,2))
  print(adf.test(dx[,1])) 
  Acf(dx[,1], main='') # ACF
  Acf(dx[,1], type="partial", main='') #PACF  
  par(mfrow=c(1,1))
}

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx[, 1]
## Dickey-Fuller = -2.5124, Lag order = 2, p-value = 0.3772
## alternative hypothesis: stationary

Interpretation:

Stationarity can be defined in precise mathematical terms, but for our purpose we mean a flat looking series, without trend, constant variance over time, a constant autocorrelation structure over time and no periodic fluctuations (seasonality). So the high p value means we do have non-stationarity. Similarly, the augmented Dickey–Fuller (ADF) statistic, used in the test, is a negative number. The more negative it is, the stronger the rejection of the null hypothesis. This is not negative enough.

Finally, the lag length (lag order) is how many terms back down the AR process you want to test for serial correlation. This is a non-stationary with AR(2).

11. Detrending and Dedrifting Data. Detrending is removing a trend from a time series; a trend usually refers to a change in the mean over time, but the overall goal is to get to stationarity. When you detrend data, you remove an aspect from the data that you think is causing some kind of distortion. For example, you might detrend data that shows an overall increase, in order to see subtrends. Usually, these subtrends are seen as fluctuations on a time series graph. To remove the linear trend, differencing is equivalent to applying a linear regression on time - “regress out” covariates.

Differencing is when a new series is constructed where the value at the current time step is calculated as the difference between the original observation and the observation at the previous time step.

detrend = function(dx) {
  dx$diff <- diff(dx[,1]) # "regressing out"
  dx <- dx[-1,]
  plot(dx$diff)
  
  par(mfrow=c(1,2))
  print(adf.test(dx$diff)) 
  Acf(dx$diff, main='') # ACF
  Acf(dx$diff, type="partial", main='') #PACF 
  par(mfrow=c(1,1))
}

dedrift = function(dx) {
  dx$logv <- log(dx[,1])
  dx <- dx[-1,]
  plot(dx$logv)
  
  par(mfrow=c(1,2))
  print(adf.test(dx$logv))
  Acf(dx$logv, main='') # ACT
  Acf(dx$logv, type="partial", main='') #PACF 
  par(mfrow=c(1,1))
}

dedriftrend = function(dx) {
  dx$logdiff <- diff(log(dx[,1]))
  dx <- dx[-1,]
  plot(dx$logdiff)
  
  par(mfrow=c(1,2))
  print(adf.test(dx$logdiff)) 
  Acf(dx$logdiff, main='') # ACT
  Acf(dx$logdiff, type="partial", main='') #PACF 
  par(mfrow=c(1,1))
}

detrend(wing_mm) # this is stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$diff
## Dickey-Fuller = -3.7455, Lag order = 2, p-value = 0.03961
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Interpretation: This time series is stationary! And as a reminder, a data that is just noise will be stationary.

12. You can also take the log() of the time series in order to stable the variance of the time series if there is a trend in the variance. Stationarity was achieved without having to stablize the variance, but you can still check:

dedrift(wing_mm) # this is stationary

## Warning in adf.test(dx$logv): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logv
## Dickey-Fuller = -6.4449, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

dedriftrend(wing_mm) # this is stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logdiff
## Dickey-Fuller = -4.2683, Lag order = 2, p-value = 0.0139
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Wing2body

# Get only bugs with long wings
data_long<-raw_data[raw_data$w_morph=="L",]

# Calculate wing2body ratio for bugs with long wings 
data_long$wing2body <- data_long$wing/as.numeric(data_long$body)

# remove NA dates
d = data_long %>%
  filter(!is.na(wing2body), !is.na(datetime))

# get wing2body averages using vectorization
ratio_avg = tapply(X=d$wing2body, INDEX=d$datetime, FUN=mean, na.rm=T)
ratio_avg = ratio_avg[!is.na(ratio_avg)]
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

ratio = xts(ratio_avg, date)
colnames(ratio) <- "wing2body"

function() {plot(ratio)}

## function() {plot(ratio)}

addEventLines(events, pos=2, srt=90, col="red")

adf.test(ratio$wing2body) # this is not stationary

## Warning in adf.test(ratio$wing2body): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ratio$wing2body
## Dickey-Fuller = 10.559, Lag order = 2, p-value = 0.99
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(ratio$wing2body, main='') # ACT
Acf(ratio$wing2body, type="partial", main='') #PACF

detrend(ratio) # this is stationary

## Warning in adf.test(dx$diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$diff
## Dickey-Fuller = -7.7839, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

dedrift(ratio) # this is not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logv
## Dickey-Fuller = -1.0499, Lag order = 1, p-value = 0.9133
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

dedriftrend(ratio)

## Warning in adf.test(dx$logdiff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logdiff
## Dickey-Fuller = -7.6591, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Wing Morph Frequency

raw_data$w_morph_binom <- NA
raw_data$wing_morph_binom[raw_data$w_morph=="S"]<-0
raw_data$wing_morph_binom[raw_data$w_morph=="L"]<-1

# remove NA dates and wing morph (S=0, L=1)
d = raw_data %>%
  filter(!is.na(wing_morph_binom), !is.na(datetime))

# get wing frequency averages using vectorization
freq_avg = tapply(X=d$wing_morph_binom, INDEX=d$datetime, FUN=mean, na.rm=T)
freq_avg = freq_avg[!is.na(freq_avg)]
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

morph = xts(freq_avg, date)
colnames(morph) <- "wing_morph_freq"

function() {plot(morph)}

## function() {plot(morph)}

addEventLines(events, pos=2, srt=90, col="red")

adf.test(morph$wing_morph_freq) # this is not stationarity

## 
##  Augmented Dickey-Fuller Test
## 
## data:  morph$wing_morph_freq
## Dickey-Fuller = -1.8526, Lag order = 2, p-value = 0.6285
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(morph$wing_morph_freq, main='') # ACT
Acf(morph$wing_morph_freq, type="partial", main='') #PACF

detrend(morph) # this is not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$diff
## Dickey-Fuller = -1.4413, Lag order = 2, p-value = 0.7852
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

dedrift(morph) # this is stationary

## Warning in adf.test(dx$logv): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logv
## Dickey-Fuller = -27.979, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

dedriftrend(morph) # this is not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logdiff
## Dickey-Fuller = -2.0047, Lag order = 2, p-value = 0.5706
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Time Series: Grouped by Sex or Host Plant

females = raw_data[raw_data$sex=="F",]
males = raw_data[raw_data$sex=="M",]

Females

Wing Length

# remove NA dates
d = females %>%
  filter(!is.na(wing), !is.na(datetime))

# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

## function() {plot(wing_mmf)}

## 
##  Augmented Dickey-Fuller Test
## 
## data:  wing_mmf$wing
## Dickey-Fuller = -2.8815, Lag order = 2, p-value = 0.2366
## alternative hypothesis: stationary

Looks like there’s not a linear trend, but changes in variation over time.

dedrift(wing_mmf)

## Warning in adf.test(dx$logv): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logv
## Dickey-Fuller = -5.4041, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Wing2body

# remove NA dates
d = data_long %>%
  filter(!is.na(wing2body), !is.na(datetime))

d = d[d$sex=="F",]

# get wing2body averages using vectorization
ratio_avg = tapply(X=d$wing2body, INDEX=d$datetime, FUN=mean, na.rm=T)
ratio_avg = ratio_avg[!is.na(ratio_avg)]
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

ratiof = xts(ratio_avg, date)
colnames(ratiof) <- "wing2body"

function() {plot(ratiof)}

## function() {plot(ratiof)}

addEventLines(events, pos=2, srt=90, col="red")

adf.test(ratiof$wing2body) # this is not stationarity

## Warning in adf.test(ratiof$wing2body): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ratiof$wing2body
## Dickey-Fuller = 1.4005, Lag order = 2, p-value = 0.99
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(ratiof$wing2body, main='') # ACT
Acf(ratiof$wing2body, type="partial", main='') #PACF

dedrift(ratiof) # still not stationary

## Warning in adf.test(dx$logv): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logv
## Dickey-Fuller = 0.051299, Lag order = 1, p-value = 0.99
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

dedriftrend(ratiof)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logdiff
## Dickey-Fuller = -2.4151, Lag order = 1, p-value = 0.4143
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red") # *still* not stationary

Wing Morph Frequency

raw_data$w_morph_binom <- NA
raw_data$wing_morph_binom[raw_data$w_morph=="S"]<-0
raw_data$wing_morph_binom[raw_data$w_morph=="L"]<-1

# remove NA dates and wing morph (S=0, L=1)
d = raw_data %>%
  filter(!is.na(wing_morph_binom), !is.na(datetime))

d = d[d$sex=="F",]
# get wing frequency averages using vectorization
freq_avg = tapply(X=d$wing_morph_binom, INDEX=d$datetime, FUN=mean, na.rm=T)
freq_avg = freq_avg[!is.na(freq_avg)]
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

morphf = xts(freq_avg, date)
colnames(morphf) <- "wing_morph_freq"

function() {plot(morphf)}

## function() {plot(morphf)}

addEventLines(events, pos=2, srt=90, col="red")

adf.test(morphf$wing_morph_freq) # this is not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  morphf$wing_morph_freq
## Dickey-Fuller = -0.77584, Lag order = 2, p-value = 0.952
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(morphf$wing_morph_freq, main='') # ACT
Acf(morphf$wing_morph_freq, type="partial", main='') #PACF

detrend(morphf) # still not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$diff
## Dickey-Fuller = -0.69068, Lag order = 2, p-value = 0.9591
## alternative hypothesis: stationary

dedriftrend(morphf) # still not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logdiff
## Dickey-Fuller = -0.83913, Lag order = 2, p-value = 0.9442
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Males

Wing Length

# remove NA dates
d = males %>%
  filter(!is.na(wing), !is.na(datetime))

# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

## function() {plot(wing_mmm)}

## 
##  Augmented Dickey-Fuller Test
## 
## data:  wing_mmm$wing
## Dickey-Fuller = -0.83612, Lag order = 2, p-value = 0.9447
## alternative hypothesis: stationary

dedrift(wing_mmm) # not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logv
## Dickey-Fuller = -1.4491, Lag order = 2, p-value = 0.7822
## alternative hypothesis: stationary

dedriftrend(wing_mm) # now stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logdiff
## Dickey-Fuller = -4.2683, Lag order = 2, p-value = 0.0139
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Wing2body

# remove NA dates
d = data_long %>%
  filter(!is.na(wing2body), !is.na(datetime))

d = d[d$sex=="M",]

# get wing2body averages using vectorization
ratio_avg = tapply(X=d$wing2body, INDEX=d$datetime, FUN=mean, na.rm=T)
ratio_avg = ratio_avg[!is.na(ratio_avg)]
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

ratiom = xts(ratio_avg, date)
colnames(ratiom) <- "wing2body"

function() {plot(ratiom)}

## function() {plot(ratiom)}

addEventLines(events, pos=2, srt=90, col="red")

adf.test(ratiom$wing2body) # this is not stationarity

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ratiom$wing2body
## Dickey-Fuller = -0.44044, Lag order = 2, p-value = 0.9776
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(ratiom$wing2body, main='') # ACT
Acf(ratiom$wing2body, type="partial", main='') #PACF

detrend(ratiom) # this is stationary

## Warning in adf.test(dx$diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$diff
## Dickey-Fuller = -7.954, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Wing Morph Frequency

raw_data$w_morph_binom <- NA
raw_data$wing_morph_binom[raw_data$w_morph=="S"]<-0
raw_data$wing_morph_binom[raw_data$w_morph=="L"]<-1

# remove NA dates and wing morph (S=0, L=1)
d = raw_data %>%
  filter(!is.na(wing_morph_binom), !is.na(datetime))

d = d[d$sex=="M",]
# get wing frequency averages using vectorization
freq_avg = tapply(X=d$wing_morph_binom, INDEX=d$datetime, FUN=mean, na.rm=T)
freq_avg = freq_avg[!is.na(freq_avg)]
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

morphm = xts(freq_avg, date)
colnames(morphm) <- "wing_morph_freq"

function() {plot(morphm)}

## function() {plot(morphm)}

addEventLines(events, pos=2, srt=90, col="red")

adf.test(morphm$wing_morph_freq) # this is not stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  morphm$wing_morph_freq
## Dickey-Fuller = -1.3425, Lag order = 2, p-value = 0.8229
## alternative hypothesis: stationary

par(mfrow=c(1,2))
Acf(morphm$wing_morph_freq, main='') # ACT
Acf(morphm$wing_morph_freq, type="partial", main='') #PACF

dedrift(morphm) # close! But marginally stationary

## 
##  Augmented Dickey-Fuller Test
## 
## data:  dx$logv
## Dickey-Fuller = -3.5464, Lag order = 2, p-value = 0.05744
## alternative hypothesis: stationary

addEventLines(events, pos=2, srt=90, col="red")

Balloon Vine

Wing Length

Wing Length

# remove NA dates
d = males %>%
  filter(!is.na(wing), !is.na(datetime))

# get wing length averages using vectorization
wing_avg = tapply(X=d$wing, INDEX=d$datetime, FUN=mean, na.rm=T)
monyear = unique(d$datetime)

# generate datetime object | datetime object needs a date, which I initialized at 01 for each month
monyeardate <- paste(monyear," 01",sep="")
date = as.Date(monyeardate, "%b %Y %d")

## function() {plot(wing_mmm)}

## 
##  Augmented Dickey-Fuller Test
## 
## data:  wing_mmm$wing
## Dickey-Fuller = -0.83612, Lag order = 2, p-value = 0.9447
## alternative hypothesis: stationary